
#include "cvf_DefaultFilters.h"

// [*] CImageFilter_Gray

CImageFilter* CImageFilter_Gray::myClone() const
{
	CImageFilter_Gray* ret = new CImageFilter_Gray(*this);
	ret->Copy_IO_Settings(this);
	return ret;
}

int CImageFilter_Gray::GetResponseDimension()
{
	return 1;
}

void CImageFilter_Gray::_ComputeResponse(int W, int H, double** ppInputResponse, double* pResponse)
{
	int j;
	int WH = W * H;
	const double* qd_R = ppInputResponse[0];
	const double* qd_G = ppInputResponse[1];
	const double* qd_B = ppInputResponse[2];
	double* pd = pResponse;
	for ( j = 0; j < WH; ++j )
	{
		*pd = 0.299 * *qd_R + 0.587 * *qd_G + 0.114 * *qd_B;
		++pd; ++qd_R; ++qd_G; ++qd_B;
	}
}

// [*] CImageFilter_Lab

CImageFilter* CImageFilter_Lab::myClone() const
{
	CImageFilter_Lab* ret = new CImageFilter_Lab(*this);
	ret->Copy_IO_Settings(this);
	return ret;
}

int CImageFilter_Lab::GetResponseDimension()
{
	return 3;
}

void CImageFilter_Lab::_ComputeResponse(int W, int H, double** ppInputResponse, double* pResponse)
{
	int j;
	int WH = W * H;
	const double* qd_R = ppInputResponse[0];
	const double* qd_G = ppInputResponse[1];
	const double* qd_B = ppInputResponse[2];
	double* pd;
	const double* qd;

	CImg<double> imgTmp;
	imgTmp.assign(W, H, 1, 3);
	pd = imgTmp.data;
	for ( j = 0; j < WH; ++j ) *(pd++) = *(qd_R++);
	for ( j = 0; j < WH; ++j ) *(pd++) = *(qd_G++);
	for ( j = 0; j < WH; ++j ) *(pd++) = *(qd_B++);
	imgTmp.RGBtoLab();

	qd = imgTmp.data;
	pd = pResponse;
	for ( j = WH*3; j > 0; --j )
		*(pd++) = *(qd++);
}

// [*] CImageFilter_Func_Sqrt

CImageFilter* CImageFilter_Func_Sqrt::myClone() const
{
	CImageFilter_Func_Sqrt* ret = new CImageFilter_Func_Sqrt(*this);
	ret->Copy_IO_Settings(this);
	return ret;
}

int CImageFilter_Func_Sqrt::GetResponseDimension()
{
	return GetInputDimension();
}

void CImageFilter_Func_Sqrt::_ComputeResponse(int W, int H, double** ppInputResponse, double* pResponse)
{
	size_t WH = (size_t)W * (size_t)H;
	int i, ii;
	size_t j;
	double* pd;
	const double* qd;

	ii = GetInputDimension();
	pd = pResponse;

	for ( i = 0; i < ii; ++i )
	{
		qd = ppInputResponse[i];
		for ( j = 0; j < WH; ++j, ++pd, ++qd )
			*pd = sqrt((*qd<EPS ? 0.0 : *qd));
	}
}

// [*] CImageFilter_Kernel2D

size_t CImageFilter_Kernel2D::GetMemoryUsage() const
{
	return CImageFilter::GetMemoryUsage() + Ker.GetMemoryUsage();
}

CImageFilter* CImageFilter_Kernel2D::myClone() const
{
	CImageFilter_Kernel2D* ret = new CImageFilter_Kernel2D(*this);
	ret->Copy_IO_Settings(this);
	return ret;
}

int CImageFilter_Kernel2D::GetResponseDimension()
{
	return GetInputDimension();
}

void CImageFilter_Kernel2D::_ComputeResponse(int W, int H, double** ppInputResponse, double* pResponse)
{
	size_t WH = (size_t)W * (size_t)H;
	int i, ii;
	double* pd;

	ii = GetInputDimension();
	pd = pResponse;

	for ( i = 0; i < ii; ++i )
	{
		Ker.Convolve(W, H, pd, ppInputResponse[i]);
		pd += WH;
	}
}

// [*] CImageFilter_Gaussian

size_t CImageFilter_Gaussian::GetMemoryUsage() const
{
	return CImageFilter::GetMemoryUsage() + Ker.GetMemoryUsage();
}

CImageFilter* CImageFilter_Gaussian::myClone() const
{
	CImageFilter_Gaussian* ret = new CImageFilter_Gaussian(*this);
	ret->Copy_IO_Settings(this);
	return ret;
}

int CImageFilter_Gaussian::GetResponseDimension()
{
	return GetInputDimension();
}

void CImageFilter_Gaussian::_ComputeResponse(int W, int H, double** ppInputResponse, double* pResponse)
{
	size_t WH = (size_t)W * (size_t)H;
	int i, ii;
	double* pd;

	ii = GetInputDimension();
	pd = pResponse;

	for ( i = 0; i < ii; ++i )
	{
		Ker.Convolve(W, H, pd, ppInputResponse[i]);
		pd += WH;
	}
}

// [*] CImageFilter_LoG

size_t CImageFilter_LoG::GetMemoryUsage() const
{
	return CImageFilter::GetMemoryUsage() + Ker.GetMemoryUsage();
}

CImageFilter* CImageFilter_LoG::myClone() const
{
	CImageFilter_LoG* ret = new CImageFilter_LoG(*this);
	ret->Copy_IO_Settings(this);
	return ret;
}

int CImageFilter_LoG::GetResponseDimension()
{
	return GetInputDimension();
}

void CImageFilter_LoG::_ComputeResponse(int W, int H, double** ppInputResponse, double* pResponse)
{
	size_t WH = (size_t)W * (size_t)H;
	int i, ii;
	double* pd;

	ii = GetInputDimension();
	pd = pResponse;

	for ( i = 0; i < ii; ++i )
	{
		Ker.Convolve(W, H, pd, ppInputResponse[i]);
		pd += WH;
	}
}

// [*] CImageFilter_1DoG

size_t CImageFilter_1DoG::GetMemoryUsage() const
{
	return CImageFilter::GetMemoryUsage() + Ker.GetMemoryUsage();
}

CImageFilter* CImageFilter_1DoG::myClone() const
{
	CImageFilter_1DoG* ret = new CImageFilter_1DoG(*this);
	ret->Copy_IO_Settings(this);
	return ret;
}

int CImageFilter_1DoG::GetResponseDimension()
{
	return GetInputDimension();
}

void CImageFilter_1DoG::_ComputeResponse(int W, int H, double** ppInputResponse, double* pResponse)
{
	size_t WH = (size_t)W * (size_t)H;
	int i, ii;
	double* pd;

	ii = GetInputDimension();
	pd = pResponse;

	for ( i = 0; i < ii; ++i )
	{
		Ker.Convolve(W, H, pd, ppInputResponse[i]);
		pd += WH;
	}
}

// [*] CImageFilter_DooG

size_t CImageFilter_DooG::GetMemoryUsage() const
{
	return CImageFilter::GetMemoryUsage() + CExternalTools::GetMemoryUsage(Ker);
}

void CImageFilter_DooG::BuildKer(CImg<double>& imgDstKer, double Sigma, double Orientation, double Aspect)
{
	double	*KerVal;
	double	center, theta, CosTheta, SinTheta, ratio, x, y;
	double	Sum;
	int		FilterSize;
	int		pos;

	FilterSize	= (int) (Sigma * 4 * Aspect) + 1;
	FilterSize	= (FilterSize % 2 == 0) ? FilterSize+1 : FilterSize;
	center		= Sigma * 2 * Aspect;
	theta		= Orientation; //(PI * Orientation)/180;    // Orientation in radians
	CosTheta	= cos(theta);
	SinTheta	= sin(theta);
	KerVal		= new double[FilterSize*FilterSize];
	ratio		= Aspect;

	for ( int i = 0; i < FilterSize; i++ )
	{
		for ( int j = 0; j < FilterSize; j++ )
		{
			x	= ((i-center) * CosTheta + (j-center) * SinTheta) / Sigma;
			y	= ((center-i) * SinTheta + (j-center) * CosTheta) / Sigma;
			pos	= FilterSize * FilterSize - 1 - (i * FilterSize + j);

			KerVal[pos]	= (double) ((CMyMath::gaussian(x - 1, 0.0, 1.0) - CMyMath::gaussian(x + 1, 0.0, 1.0)) * CMyMath::gaussian(y, 0.0, ratio));
		}
	}

	Sum	= 0;
	for ( int i = 0; i < FilterSize * FilterSize; i++ )
		Sum	+= KerVal[i];
	Sum	/= FilterSize * FilterSize;
	for ( int i = 0; i < FilterSize * FilterSize; i++ )
		KerVal[i]	-= Sum;
	Sum	= 0;
	for ( int i = 0; i < FilterSize * FilterSize; i++ )
		Sum	+= fabs(KerVal[i]);
	for ( int i = 0; i < FilterSize * FilterSize; i++ )
		KerVal[i]	/= Sum;

	imgDstKer.assign(FilterSize, FilterSize);
	for ( int i = 0; i < FilterSize; i++ )
	{
		for ( int j = 0; j < FilterSize; j++ )
		{
			pos = i * FilterSize + j;
			imgDstKer.at(i, j)	= KerVal[pos];
		}
	}

	delete[] KerVal;
}

CImageFilter* CImageFilter_DooG::myClone() const
{
	CImageFilter_DooG* ret = new CImageFilter_DooG(*this);
	ret->Copy_IO_Settings(this);
	return ret;
}

int CImageFilter_DooG::GetResponseDimension()
{
	return GetInputDimension();
}

void CImageFilter_DooG::_ComputeResponse(int W, int H, double** ppInputResponse, double* pResponse)
{
	int WH = W * H;
	int i, j, ii;
	const double* qd;
	double* pd;
	CImg<double> tImgSrc;
	CImg<double> tImgDst;
	CImg<int> imgOutput;

	int _j;
	CImg<double> tImgTmpA;
	CImg<double> tImgTmpB;
	CTransformableImage tfmImage_a;
	CImageTransformation a_tfm;
	CImageTransformation b_tfm;
	double nw, nh;
	int iw, ih, iwh;
	double tx, ty;
	CKernel2D<double> tKer;

	tImgDst.assign(W, H, 1, 1);
	CKernel2D<double>::GenerateGaussian(tKer, Sigma, Sigma * Aspect);

	ii = GetInputDimension();
	pd = pResponse;
	for ( i = 0; i < ii; ++i )
	{
		tImgSrc.assign(ppInputResponse[i], W, H, 1, 1, true);

		if (1)
		{
			tImgDst.fill(0.0);

			//CTimer ta;
			//ta.Reset();
			//ta.Start();

			tfmImage_a.myInit(tImgSrc);
			tfmImage_a.DefaultColor.Fill(128.0);
			tfmImage_a.OutExtent = Sigma * 4.0 + 0.5;
			a_tfm.Clear();
			a_tfm.PushBack(new C_ImgTran_Rotate(W * 0.5, H * 0.5, Orientation));
			tfmImage_a.NormalizeTransformedSize(nw, nh, a_tfm);
			iw = (int)(nw + 0.5); ih = (int)(nh + 0.5); iwh = iw * ih;
			b_tfm = a_tfm;
			b_tfm.ToInverseTransformation();

			//ta.Stop(); cout<<"d "<<setprecision(2)<<fixed<<ta.GetElapsedTime()<<endl;
			//ta.Start();

			tImgTmpA.assign(iw, ih, 1, 1);
			for ( _j = 0; _j < iwh; ++_j )
			{
				a_tfm.Inverse(tx, ty, _j%iw + 0.5 + Sigma, _j/iw + 0.5);
				tImgTmpA.data[_j] = tfmImage_a.GetPixel(tx, ty, 0, CDenseImage::EXPANDING_MODE_EDGE);

				a_tfm.Inverse(tx, ty, _j%iw + 0.5 - Sigma, _j/iw + 0.5);
				tImgTmpA.data[_j] -= tfmImage_a.GetPixel(tx, ty, 0, CDenseImage::EXPANDING_MODE_EDGE);
			}

			//ta.Stop(); cout<<"d "<<setprecision(2)<<fixed<<ta.GetElapsedTime()<<endl;
			//ta.Start();

			tImgTmpB.assign(iw, ih, 1, 1);
			tKer.Convolve(iw, ih, tImgTmpB.data, tImgTmpA.data);

			//ta.Stop(); cout<<"d "<<setprecision(2)<<fixed<<ta.GetElapsedTime()<<endl;
			//ta.Start();

			tfmImage_a.myInit(tImgTmpB);
			tfmImage_a.BuildImage(tImgDst.data, W, H, b_tfm, CDenseImage::EXPANDING_MODE_COLOR);

			//ta.Stop(); cout<<"d "<<setprecision(2)<<fixed<<ta.GetElapsedTime()<<endl;
		}
		else
			tImgDst = tImgSrc.get_convolve(Ker);

		//CExternalTools::NormalizeToGeneralIntGrayscaleImage(tImgDst, imgOutput); imgOutput.display();

		qd = tImgDst.data;
		for ( j = 0; j < WH; ++j )
			*(pd++) = *(qd++);
	}
}

// [*] CImageFilter_DoG

size_t CImageFilter_DoG::GetMemoryUsage() const
{
	return CImageFilter::GetMemoryUsage() + CExternalTools::GetMemoryUsage(Ker);
}

void CImageFilter_DoG::BuildKer(CImg<double>& imgDstKer, double Sigma)
{
	double	*KerVal;
	double	Sigma1, Sigma2, Sigma3;
	int		center;
	double	x, y;
	double	Sum;
	int		FilterSize;
	int		pos;

	Sigma1		= 0.62;
	Sigma2		= 1.0;
	Sigma3		= 1.6;
	FilterSize	= (int) (Sigma * 6);
	FilterSize	= (FilterSize % 2 == 0) ? FilterSize+1 : FilterSize;
	center		= int(FilterSize / 2 + 0.5);
	KerVal		= new double[FilterSize*FilterSize];

	//?!

	for ( int i = 0; i < FilterSize; i++ )
	{
		for ( int j = 0; j < FilterSize; j++ )
		{
			x	= (i - center) / Sigma;
			y	= (j - center) / Sigma;
			pos	= FilterSize * FilterSize - 1 - (i * FilterSize + j);

			KerVal[pos]	= (double) (- (CMyMath::gaussian(x, 0.0, Sigma1) * CMyMath::gaussian(y, 0.0, Sigma1)));
			KerVal[pos]	= (double) (- (CMyMath::gaussian(x, 0.0, Sigma3) * CMyMath::gaussian(y, 0.0, Sigma3))) + KerVal[pos];
			KerVal[pos]	= (double) (2 * (CMyMath::gaussian(x, 0.0, Sigma2) * CMyMath::gaussian(y, 0.0, Sigma2))) + KerVal[pos];
		}
	}

	Sum	= 0;
	for ( int i = 0; i < FilterSize * FilterSize; i++ )
		Sum	+= KerVal[i];
	Sum	/= FilterSize * FilterSize;
	for ( int i = 0; i < FilterSize * FilterSize; i++ )
		KerVal[i]	-= Sum;
	Sum	= 0;
	for ( int i = 0; i < FilterSize * FilterSize; i++ )
		Sum	+= fabs(KerVal[i]);
	for ( int i = 0; i < FilterSize * FilterSize; i++ )
		KerVal[i]	/= Sum;

	imgDstKer.assign(FilterSize, FilterSize);
	for ( int i = 0; i < FilterSize; i++ )
	{
		for ( int j = 0; j < FilterSize; j++ )
		{
			pos				= i * FilterSize + j;
			imgDstKer.at(i, j)	= KerVal[pos];
		}
	}

	delete[] KerVal;
}

CImageFilter* CImageFilter_DoG::myClone() const
{
	CImageFilter_DoG* ret = new CImageFilter_DoG(*this);
	ret->Copy_IO_Settings(this);
	return ret;
}

int CImageFilter_DoG::GetResponseDimension()
{
	return GetInputDimension();
}

void CImageFilter_DoG::_ComputeResponse(int W, int H, double** ppInputResponse, double* pResponse)
{
	int WH = W * H;
	int i, j, ii;
	const double* qd;
	double* pd;
	CImg<double> tImgSrc;
	CImg<double> tImgDst;

	ii = GetInputDimension();
	pd = pResponse;
	for ( i = 0; i < ii; ++i )
	{
		tImgSrc.assign(ppInputResponse[i], W, H, 1, 1, true);
		tImgDst = tImgSrc.get_convolve(Ker);
		qd = tImgDst.data;
		for ( j = 0; j < WH; ++j )
			*(pd++) = *(qd++);
	}
}

// [*] CImageFilter_Gradient

size_t CImageFilter_Gradient::GetMemoryUsage() const
{
	return CImageFilter::GetMemoryUsage() + sizeof(double) + sizeof(bool);
}

CImageFilter_Gradient::CImageFilter_Gradient(double GaussianSigma, bool fUseSignedOrientation)
{
	this->GaussianSigma = GaussianSigma;
	this->fUseSignedOrientation = fUseSignedOrientation;
}

CImageFilter* CImageFilter_Gradient::myClone() const
{
	CImageFilter_Gradient* ret = new CImageFilter_Gradient(*this);
	ret->Copy_IO_Settings(this);
	return ret;
}

int CImageFilter_Gradient::GetResponseDimension()
{
	return 4;
}

void CImageFilter_Gradient::_ComputeResponse(int W, int H, double** ppInputResponse, double* pResponse)
{
	size_t WH = ((size_t)W) * (size_t)H;
	int i,j,ii;
	double* pd;
	double* p_mag;
	double* p_ori;
	double* p_gx;
	double* p_gy;
	CSimpleTypeArray<double> max_mag;
	const double* q_mag;
	const double* q_ori;
	const double* q_gx;
	const double* q_gy;

	CImageGradient<double> ig;
	CImg<double> imgTmp;
	imgTmp.assign(W, H);

	max_mag.Resize(WH);

	ii = GetInputDimension();
	for ( i = 0; i < ii; ++i )
	{
		pd = imgTmp.data;
		q_mag = ppInputResponse[i];
		for ( j = 0; j < (int)WH; ++j )
			*(pd++) = *(q_mag++);

		if (GaussianSigma > 0.05)
		{
			CKernel2D<double> ker;
			CKernel2D<double>::GenerateGaussian(ker, GaussianSigma);
			CImg<double> imgTT;
			imgTT.assign(imgTmp, false);
			ker.Convolve(imgTT.dimx(), imgTT.dimy(), imgTmp.data, imgTT.data);
		}

		ig.InitGradient(imgTmp, CImageGradient<double>::USING_MASK_XY1x3, fUseSignedOrientation);

		pd = max_mag.pElement;
		p_mag = pResponse + WH * 0;
		p_ori = pResponse + WH * 1;
		p_gx = pResponse + WH * 2;
		p_gy = pResponse + WH * 3;

		q_mag = ig.MagnitudeImage.data;
		q_ori = ig.OrientationImage.data;
		q_gx = ig.xGradientImage.data;
		q_gy = ig.yGradientImage.data;
		for ( j = 0; j < (int)WH; ++j, ++pd, ++q_mag, ++q_ori, ++q_gx, ++q_gy )
		{
			if (i==0 || *pd<*q_mag)
			{
				*pd = *q_mag;
				*p_mag = *q_mag;
				*p_ori = *q_ori;
				*p_gx = *q_gx;
				*p_gy = *q_gy;
			}
			++p_mag; ++p_ori; ++p_gx; ++p_gy;
		}
	}
}

// [*] CImageFilter_2D_Histogram

CImageFilter_2D_Histogram::CImageFilter_2D_Histogram(
	int ResponseInputIndex_Weight, int ResponseInputIndex_Orientation,
	int BinAmount, double ValueBoundA, double ValueBoundB,
	int BinSelectMode, bool fAreBinWrapped)
{
	this->ResponseInputIndex_Weight = ResponseInputIndex_Weight;
	this->ResponseInputIndex_Orientation = ResponseInputIndex_Orientation;
	if (BinAmount<1) BinAmount = 1;
	this->BinAmount = BinAmount;
	this->ValueBoundA = ValueBoundA;
	this->ValueBoundB = ValueBoundB;
	this->BinSelectMode = BinSelectMode;
	this->fAreBinWrapped = fAreBinWrapped;
	HistogramConfig.myInit(BinAmount, fAreBinWrapped, ValueBoundA, ValueBoundB);
}

size_t CImageFilter_2D_Histogram::GetMemoryUsage() const
{
	return CImageFilter::GetMemoryUsage()
		+ sizeof(BinAmount) + sizeof(ValueBoundA) + sizeof(ValueBoundB) + sizeof(HistogramConfig)
		+ sizeof(BinSelectMode) + sizeof(fAreBinWrapped);
}

CImageFilter* CImageFilter_2D_Histogram::myClone() const
{
	CImageFilter_2D_Histogram* ret = new CImageFilter_2D_Histogram(*this);
	ret->Copy_IO_Settings(this);
	return ret;
}

int CImageFilter_2D_Histogram::GetResponseDimension()
{
	return BinAmount;
}

void CImageFilter_2D_Histogram::_ComputeResponse(int W, int H, double** ppInputResponse, double* pResponse)
{
	int ii = GetInputDimension();
	size_t WH = ((size_t)W) * (size_t)H;
	size_t WHV = WH * BinAmount;

	CSimpleTypeArray<double> FakeWeightList;
	CSimpleTypeArray<double> FakeOrientationList;
	const double* qw;
	const double* qo;

	if (ResponseInputIndex_Weight<0 || ResponseInputIndex_Weight>=ii)
	{
		ResponseInputIndex_Weight = -1;
		FakeWeightList.Resize(WH);
		FakeWeightList.Fill(1.0);
		qw = FakeWeightList.pElement;
	}
	else
		qw = ppInputResponse[ResponseInputIndex_Weight];

	if (ResponseInputIndex_Orientation<0 || ResponseInputIndex_Orientation>=ii)
	{
		ResponseInputIndex_Orientation = -1;
		FakeOrientationList.Resize(WH);
		FakeOrientationList.Fill((ValueBoundA + ValueBoundB) * 0.5);
		qo = FakeOrientationList.pElement;
	}
	else
		qo = ppInputResponse[ResponseInputIndex_Orientation];

	double* pd;
	size_t i;
	int b1, b2;
	double v;

	pd = pResponse;
	for ( i = 0; i < WHV; ++i, ++pd )
		*pd = 0.0;

	v = 0.0;
	b1 = b2 = 0;
	pd = pResponse;
	for ( i = 0; i < WH; ++i, ++pd, ++qw, ++qo )
	{
		if (BinSelectMode == MODE_BOUNDING)
		{
			v = HistogramConfig.SelectBin_Bounding(b1, b2, *qo);
			pd[WH * (size_t)b1] += (*qw) * (1-v);
			pd[WH * (size_t)b2] += (*qw) * v;
		}
		else // if (BinSelectMode == MODE_NEAREST)
		{
			b1 = HistogramConfig.SelectBin_Nearest(*qo);
			pd[WH * (size_t)b1] += *qw;
		}
	}
}

// [*] CImageFilter_LinearCombination

CImageFilter_LinearCombination::CImageFilter_LinearCombination()
{
	Clear_IO_Relationship();
}

void CImageFilter_LinearCombination::Clear_IO_Relationship()
{
	CImageFilter::Clear_IO_Relationship();
	InputWeightList.Clear();
	InputOffset_X_List.Clear();
	InputOffset_Y_List.Clear();
}

void CImageFilter_LinearCombination::PushBackInputResponse(
	CImageFilter* pFilter, int ResponseIndex,
	double Weight, int Offset_X, int Offset_Y)
{
	CImageFilter::PushBackInputResponse(pFilter, ResponseIndex);
	InputWeightList.PushBack(Weight);
	InputOffset_X_List.PushBack(Offset_X);
	InputOffset_Y_List.PushBack(Offset_Y);
}

size_t CImageFilter_LinearCombination::GetMemoryUsage() const
{
	return CImageFilter::GetMemoryUsage()
		+ InputWeightList.GetMemoryUsage()
		+ InputOffset_X_List.GetMemoryUsage()
		+ InputOffset_Y_List.GetMemoryUsage();
}

CImageFilter* CImageFilter_LinearCombination::myClone() const
{
	CImageFilter_LinearCombination* ret = new CImageFilter_LinearCombination(*this);
	ret->Copy_IO_Settings(this);
	return ret;
}

int CImageFilter_LinearCombination::GetResponseDimension()
{
	return 1;
}

void CImageFilter_LinearCombination::_ComputeResponse(int W, int H, double** ppInputResponse, double* pResponse)
{
	int WH = W * H;
	int ii = GetInputDimension();

	double Weight;
	int Offset_X, Offset_Y;
	int i, ix, iy;
	int lx, hx, ly, hy;
	double* pd;
	const double* qd;

	for ( i = 0; i < ii; ++i )
	{
		Weight = InputWeightList.pElement[i];
		Offset_X = InputOffset_X_List.pElement[i];
		Offset_Y = InputOffset_Y_List.pElement[i];

		lx = 0 - (Offset_X<0 ? Offset_X : 0);
		ly = 0 - (Offset_Y<0 ? Offset_Y : 0);
		hx = W - (Offset_X>0 ? Offset_X : 0);
		hy = H - (Offset_Y>0 ? Offset_Y : 0);

		for ( iy = ly; iy < hy; ++iy )
		{
			pd = pResponse + (lx + iy * W);
			qd = ppInputResponse[i] + ((lx + Offset_X) + (iy + Offset_Y) * W);
			for ( ix = lx; ix < hx; ++ix )
				*(pd++) += *(qd++) * Weight;
		}
	}
}

// [*] CImageFilter_CannyEdge

size_t CImageFilter_CannyEdge::GetMemoryUsage() const
{
	return CImageFilter::GetMemoryUsage() + sizeof(double) * 5 + sizeof(bool);
}

CImageFilter_CannyEdge::CImageFilter_CannyEdge(double Sigma, double tLow, double tHigh, double PrimaryOrientation, double OrientationWidth, bool fUseSignedOrientation)
{
	this->Sigma = Sigma;
	this->tLow = tLow;
	this->tHigh = tHigh;
	this->OriMin = PrimaryOrientation - OrientationWidth / 2.0;
	this->OriMax = PrimaryOrientation + OrientationWidth / 2.0;
	this->fUseSignedOrientation = fUseSignedOrientation;
}

CImageFilter* CImageFilter_CannyEdge::myClone() const
{
	CImageFilter_CannyEdge* ret = new CImageFilter_CannyEdge(*this);
	ret->Copy_IO_Settings(this);
	return ret;
}

int CImageFilter_CannyEdge::GetResponseDimension()
{
	return 2;
}

void CImageFilter_CannyEdge::_ComputeResponse(int W, int H, double** ppInputResponse, double* pResponse)
{
	int iv, va;
	int WH = W * H;

	CImg<double> tImgSrc;
	CImg<unsigned char> iImgSrc;

	CImg<int> imgDeltaX;
	CImg<int> imgDeltaY;
	CImg<float> imgOri;
	CImg<int> imgMag;

	CImg<int> _imgDeltaX;
	CImg<int> _imgDeltaY;
	CImg<float> _imgOri;
	CImg<int> _imgMag;

	CImg<unsigned char> imgEdge;

	int i;
	double t, tori;
	const double* src;
	double* pMag;
	double* pOri;
	const unsigned char* qe;
	const float* qo;

	va = GetInputDimension();

	pMag = pResponse;
	pOri = pResponse + WH;
	for ( i = 0; i < WH; ++i )
	{
		*(pMag++) = 0.0;
		*(pOri++) = 0.0;
	}

	imgEdge.assign(W, H);
	imgDeltaX.assign(W, H); imgDeltaY.assign(W, H); imgMag.assign(W, H); imgOri.assign(W, H);
	_imgDeltaX.assign(W, H); _imgDeltaY.assign(W, H); _imgMag.assign(W, H); _imgOri.assign(W, H);
	for ( iv = 0; iv < va; ++iv )
	{
		src = ppInputResponse[iv];
		tImgSrc.assign(src, W, H, 1, 1, true);
		CExternalTools::NormalizeToGeneralIntGrayscaleImage(tImgSrc, iImgSrc);

		sfu_canny_edge_detector::sfu_canny_edge_Gradient(iImgSrc.data, H, W, (float)Sigma, _imgDeltaX.data, _imgDeltaY.data, _imgOri.data, _imgMag.data);

		for ( i = 0; i < WH; ++i )
		{
			if (iv>0 && imgMag.data[i]>_imgMag.data[i]) continue;
			imgOri.data[i] = _imgOri.data[i];
			imgMag.data[i] = _imgMag.data[i];
			imgDeltaX.data[i] = _imgDeltaX.data[i];
			imgDeltaY.data[i] = _imgDeltaY.data[i];
		}
	}

	//char strTemp[10001];
	//sprintf_s(strTemp, "%d %d %.2lf %.2lf", H, W, tLow, tHigh);
	//cout<<strTemp<<endl;
	//sprintf_s(strTemp, "%d %d %d %d %d",
	//	(int)imgEdge.size(), (int)imgDeltaX.size(), (int)imgDeltaY.size(), (int)imgOri.size(), (int)imgMag.size());
	//cout<<strTemp<<endl;
	sfu_canny_edge_detector::sfu_canny_edge_NMS(H, W, (float)tLow, (float)tHigh, imgEdge.data, imgDeltaX.data, imgDeltaY.data, imgOri.data, imgMag.data);
	//cout<<"ok b"<<endl;

	//CImg<unsigned char> imgAAA;
	//CImg<int> imgOutput;
	//imgAAA.assign(imgEdge.data, W, H, 1, 1, true);
	//CExternalTools::NormalizeToGeneralIntGrayscaleImage(imgAAA, imgOutput);
	//imgOutput.display();

	pMag = pResponse;
	pOri = pResponse + WH;
	qe = imgEdge.data;
	qo = imgOri.data;
	for ( i = 0; i < WH; ++i, ++pMag, ++pOri, ++qe, ++qo )
	{
		//*pMag = 0.0;
		if (*qe != 0) continue;
		tori = 0.0 - (*qo);

		if (fUseSignedOrientation)
		{
			tori = CMyMath::FormatAngle(tori);
			t = CMyMath::FormatAngle(tori - OriMin, PI);
		}
		else
		{
			tori = CMyMath::FormatAngle_Small(tori);
			t = CMyMath::FormatAngle_Small(tori - OriMin, PI/2);
		}
		if (t < OriMax - OriMin)
		{
			*pMag = 1.0;
			*pOri = tori;
		}
	}

	//for ( iv = 0; iv < va; ++iv )
	//{
	//	src = ppInputResponse[iv];
	//	tImgSrc.assign(src, W, H, 1, 1, true);
	//	CExternalTools::NormalizeToGeneralIntGrayscaleImage(tImgSrc, iImgSrc);
	//	imgEdge.assign(W, H);
	//	imgDeltaX.assign(W, H); imgDeltaY.assign(W, H); imgMag.assign(W, H); imgOri.assign(W, H);
	//	_imgDeltaX.assign(W, H); _imgDeltaY.assign(W, H); _imgMag.assign(W, H); _imgOri.assign(W, H);

	//	sfu_canny_edge_detector::sfu_canny_image(iImgSrc.data, H, W, (float)Sigma, (float)tLow, (float)tHigh, imgEdge.data, imgOri.data, imgMag.data);

	//	pMag = pResponse;
	//	pOri = pResponse + WH;
	//	qe = imgEdge.data;
	//	qo = imgOri.data;
	//	for ( i = 0; i < WH; ++i, ++pMag, ++pOri, ++qe, ++qo )
	//	{
	//		//*pMag = 0.0;
	//		if (*qe != 0) continue;
	//		tori = CMyMath::FormatAngle(-(*qo));

	//		if (!fUseSignedOrientation)
	//			tori = CMyMath::FormatAngle_Small(tori);

	//		t = OriMin - tori;
	//		if (t > 0)
	//		{
	//			t = floor((t/2.0+PI-EPS)/PI);
	//			tori += t*(PI+PI);
	//		}
	//		t = tori - OriMax;
	//		if (t > 0)
	//		{
	//			t = floor((t/2.0+PI-EPS)/PI);
	//			tori -= t*(PI+PI);
	//		}
	//		
	//		if (tori < OriMax)
	//		{
	//			*pMag = 1.0;
	//			*pOri = tori;
	//		}
	//	}
	//}
}

// [*] CEdgeForest

int CEdgeForest::NMOV = 8;
const int CEdgeForest::XMOV[MAX_NMOV] = {-1,0,0,+1,+1,-1,+1,-1};
const int CEdgeForest::YMOV[MAX_NMOV] = {0,-1,+1,0,+1,-1,-1,+1};

void CEdgeForest::myRelease()
{
	W = H = WH = 0;
	CDataStructureTools::MyRelease_List(_matMag);
	CDataStructureTools::MyRelease_List(_matOri);
	CDataStructureTools::MyRelease_List(rgNode);
}

CEdgeForest::CEdgeForest()
{
	W = H = WH = 0;
	_matMag = NULL;
	_matOri = NULL;
	rgNode = NULL;
	myRelease();
}

CEdgeForest::~CEdgeForest()
{
	myRelease();
}

void CEdgeForest::GenerateForest(double DirectionThreshold)
{
	int i;
	CNodeOfEdgeTree* pCurrentNode;

	for ( i = 0; i < WH; ++i )
	{
		pCurrentNode = rgNode + i;

		if (_matMag[i])
		{
			pCurrentNode->Label = -1;
			pCurrentNode->OriginalDirection = _matOri[i];
		}
		else
		{
			pCurrentNode->Label = -2;
			pCurrentNode->OriginalDirection = 0.0f;
		}

		pCurrentNode->SmoothedDirection = pCurrentNode->OriginalDirection;

		pCurrentNode->pParent = NULL;
		pCurrentNode->maxEdgeLength = 0;
		pCurrentNode->maxUpLength = 0;
		pCurrentNode->maxDownLength = 0;
		pCurrentNode->pDownLength1st = NULL;
		pCurrentNode->pDownLength2nd = NULL;
	}

	_DirectionThreshold = DirectionThreshold + EPS;

	TreeCount = 0;
	for ( i = 0; i < WH; ++i )
	{
		pCurrentNode = rgNode + i;
		if (pCurrentNode->Label == -1)
		{
			TraceFirstRound(pCurrentNode);
			++TreeCount;
		}
	}

	TreeCount = 0;
	for ( i = 0; i < WH; ++i )
	{
		pCurrentNode = rgNode + i;
		if (pCurrentNode->Label>=TreeCount)
		{
			pCurrentNode->maxUpLength = 0;
			TraceSecondRound(pCurrentNode);
			++TreeCount;
		}
	}
}

void CEdgeForest::TraceFirstRound(CNodeOfEdgeTree* pRoot)
{
	int i,dl;
	int tx,ty;
	CNodeOfEdgeTree *pNode;

	pRoot->Label = TreeCount;
	pRoot->pDownLength1st = pRoot->pDownLength2nd = NULL;
	dl = 0;
	for ( i = 0; i < NMOV; ++i )
	{
		tx = pRoot->x + XMOV[i];
		ty = pRoot->y + YMOV[i];
		if (tx<0 || ty<0 || tx>=W || ty>=H) continue;
		pNode = pRoot + OffsetMov[i];

		if (pNode->Label != -1) continue;
		if (pNode == pRoot->pParent) continue;
		
		if (fabs(CMyMath::FormatAngle(pRoot->SmoothedDirection - pNode->OriginalDirection))>_DirectionThreshold) continue;
		//if (fabs(CMyMath::FormatAngle_Small(pRoot->SmoothedDirection - pNode->OriginalDirection))>_DirectionThreshold) continue;

		pNode->SmoothedDirection = CMyMath::FormatAngle(pNode->OriginalDirection + 0.5 * CMyMath::FormatAngle(pRoot->OriginalDirection - pNode->OriginalDirection));
		if (pRoot->pParent != NULL)
		{
			pNode->SmoothedDirection = CMyMath::FormatAngle(pNode->SmoothedDirection + 1.0 / 3.0 * CMyMath::FormatAngle(pRoot->pParent->OriginalDirection - pNode->SmoothedDirection));
		}

		pNode->pParent = pRoot;
		TraceFirstRound(pNode);
		if (pNode->maxDownLength > dl)
		{
			dl = pNode->maxDownLength;
			pRoot->pDownLength2nd = pRoot->pDownLength1st;
			pRoot->pDownLength1st = pNode;
		}
		else
			if ((pRoot->pDownLength2nd == NULL) || (pRoot->pDownLength2nd->maxDownLength < pNode->maxDownLength))
				pRoot->pDownLength2nd = pNode;
	}
	pRoot->maxDownLength = ++dl;
}

void CEdgeForest::TraceSecondRound(CNodeOfEdgeTree* pRoot)
{
	int i,t;
	int tx,ty;
	CNodeOfEdgeTree *pNode;

	if (pRoot->pDownLength2nd != NULL)
		t = pRoot->pDownLength2nd->maxDownLength;
	else
		t = 0;
	if (t < pRoot->maxUpLength)
		t = pRoot->maxUpLength;

	pRoot->maxEdgeLength = pRoot->maxDownLength + t;

	for ( i = 0; i < NMOV; ++i )
	{
		tx = pRoot->x + XMOV[i];
		ty = pRoot->y + YMOV[i];
		if (tx<0 || ty<0 || tx>=W || ty>=H) continue;
		pNode = pRoot + OffsetMov[i];
		if (pNode->pParent != pRoot) continue;

		if (pNode != pRoot->pDownLength1st)
			t = pRoot->pDownLength1st->maxDownLength;
		else if (pRoot->pDownLength2nd != NULL)
			t = pRoot->pDownLength2nd->maxDownLength;
		else
			t = 0;
		if (t < pRoot->maxUpLength)
			t = pRoot->maxUpLength;

		pNode->maxUpLength = ++t;

		TraceSecondRound(pNode);
	}
}

// [*] CImageFilter_EdgeForest

size_t CImageFilter_EdgeForest::GetMemoryUsage() const
{
	return CImageFilter::GetMemoryUsage() + sizeof(double) * 2 + sizeof(bool) + sizeof(int);
}

CImageFilter_EdgeForest::CImageFilter_EdgeForest(double MagThreshold, bool fTrueMagSign, int LengthThreshold, double DirectionThreshold)
{
	this->MagThreshold = MagThreshold;
	this->fTrueMagSign = fTrueMagSign;
	this->LengthThreshold = LengthThreshold;
	this->DirectionThreshold = DirectionThreshold;
}

CImageFilter* CImageFilter_EdgeForest::myClone() const
{
	CImageFilter_EdgeForest* ret = new CImageFilter_EdgeForest(*this);
	ret->Copy_IO_Settings(this);
	return ret;
}

int CImageFilter_EdgeForest::GetResponseDimension()
{
	return 2;
}

void CImageFilter_EdgeForest::_ComputeResponse(int W, int H, double** ppInputResponse, double* pResponse)
{
	int WH = W * H;

	if (1)
	{
		CEdgeForest ef;
		int i;
		double* pMag;
		double* pOri;
		const CEdgeForest::CNodeOfEdgeTree* qNode;

		ef.myInit(W, H, ppInputResponse[0], ppInputResponse[1], MagThreshold, fTrueMagSign);
		ef.GenerateForest(DirectionThreshold);

		pMag = pResponse; pOri = pResponse + WH;
		qNode = ef.rgNode;
		for ( i = 0; i < WH; ++i, ++pMag, ++pOri, ++qNode )
		{
			if (qNode->maxEdgeLength < LengthThreshold)
			{
				*pMag = 0.0;
				*pOri = 0.0;
			}
			else
			{
				*pMag = 1.0;
				*pOri = qNode->OriginalDirection;
			}
		}

		//CImg<double> imgAAA;
		//CImg<int> imgOutput;
		//imgAAA.assign(pResponse, W, H, 1, 1, true);
		//CExternalTools::NormalizeToGeneralIntGrayscaleImage(imgAAA, imgOutput);
		//imgOutput.display();
	}
}

////////////////////////////////////////////////////////////////

